!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!! Error computations.
!!
!! \n
!!
!! All the subroutines of this module assume that an analytic solution
!! is available. Errors are computed with numerical quadrature, and
!! it's advisable to use more accurate quadrature formulas than those
!! used in finite element method. The extra accuracy is controlled by
!! the input argument \c extra_deg used by all the error computation
!! subroutines.
!<----------------------------------------------------------------------
module mod_error_norms

!-----------------------------------------------------------------------

 !$ use omp_lib

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_sympoly, only: &
   mod_sympoly_initialized, &
   t_sympol,     &
   operator(*),  &
   operator(**), &
   me_int,       &
   ev

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid, affmap, &
   diameter

 use mod_base, only: &
   mod_base_initialized, &
   t_base, &
   new_base => base

 use mod_testcases, only: &
   mod_testcases_initialized, &
   i_coeff_lam, i_coeff_xiadv,&
   i_coeff_diff, i_coeff_adv, &
   i_coeff_q, i_coeff_divq

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_error_norms_constructor, &
   mod_error_norms_destructor,  &
   mod_error_norms_initialized, &
   hybrid_err, &
   primal_err, &
   primal_err_sef, &
   dual_err,   &
   primal_flux_err

 private

!-----------------------------------------------------------------------

! Module types and parameters

! Module variables

 ! public members
 logical, protected ::               &
   mod_error_norms_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_error_norms'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_error_norms_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_sympoly_initialized.eqv..false.)  .or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_base_initialized.eqv..false.)     .or. &
       (mod_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_error_norms_initialized.eqv..true.) then
   endif
   !----------------------------------------------

   mod_error_norms_initialized = .true.
 end subroutine mod_error_norms_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_error_norms_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_error_norms_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_error_norms_initialized = .false.
 end subroutine mod_error_norms_destructor

!-----------------------------------------------------------------------
 
 !> Error computation for the hybrid variable.
 !!
 !! The following errors are computed:
 !! <ul>
 !!  <li>\c e_l2s
 !! \f$=\left( \sum_s |s|\|u - u_h\|^2 \right)^{\frac{1}{2}} \f$
 !!  <li>\c e_l2k
 !! \f$=\left(\sum_K\sum_{s_K}h_K\|u - u_h\|^2\right)^{\frac{1}{2}}\f$
 !!  <li>\c e_pl2s
 !! \f$=\left( \sum_s |s|\|\Pi_hu - u_h\|^2 \right)^{\frac{1}{2}} \f$
 !!  <li>\c e_pl2k
 !! \f$=\left(\sum_K\sum_{s_K}h_K\|\Pi_hu-u_h\|^2\right)^{\frac{1}{2}}\f$
 !! </ul>
 !! where \f$s\f$ and \f$K\f$ denote sides and elements, respectively,
 !! \f$\|u\|^2 = \int_s u^2\,d\sigma\f$, \f$h_K = \max_{s \in \partial
 !! K}{ |s| }\f$ and \f$\Pi_h\f$ is the \f$L^2\f$ projector (computed
 !! exploiting the orthogonality of the side basis) 
 !<
 subroutine hybrid_err(grid,base_in,lamh,lam,extra_deg, &
                       e_l2s,e_l2k,e_pl2s,e_pl2k)
  type(t_grid), target, intent(in) :: grid !< grid
  type(t_base), intent(in) :: base_in !< original basis
  real(wp), intent(in) :: lamh(:)     !< numerical solution
  procedure(i_coeff_lam) :: lam       !< exact solution
  integer, intent(in) :: extra_deg    !< extra precision
  real(wp), intent(out) :: e_l2s, e_l2k, e_pl2s, e_pl2k
 
  integer :: i, is
  real(wp) :: l2s, pl2s
  real(wp) ::  hs, hk1, hk2
  real(wp), allocatable :: e_norm(:,:), xgs(:,:), wgs(:), lams(:), &
    plams(:), lamhs(:), lamk(:), lamhk(:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=200) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'hybrid_err'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , e_s=base_in%e_s , &
                    degs=base_in%degs+extra_deg )

   ! normalization (even though usually base is already normalized)
   allocate( e_norm(base%nk,base%ms) )
   do i=1,base%nk
     hs = sqrt(me_int(base%me%d-1,base%e_s(i)**2)) ! hs used as temporary
     e_norm(i,:) = ev( (1.0_wp/hs)*base%e_s(i) , base%xigs )
   enddo

   ! omp directives commented because of a bug in ifort
   ! this is probably the bug "ifort-omp-ice-with-interface"
   ! !$omp parallel private(xgs,wgs,lams,plams,lamhs,lamk,lamhk,&
   ! !$omp                  is,as,hk1,hk2,i, l2s,pl2s)            &
   ! !$omp           shared(base,grid,lamh,e_norm,ndim,           &
   ! !$omp                  e_l2s,e_l2k,e_pl2s,e_pl2k)
   ! !$omp          default(none)

   allocate( xgs(grid%m,base%ms), wgs(base%ms),  &
     lams(base%ms), plams(base%ms), lamhs(base%ms), &
     lamk(base%nk), lamhk(base%nk) )

   ! !$omp single
    e_l2s  = 0.0_wp
    e_l2k  = 0.0_wp
    e_pl2s = 0.0_wp
    e_pl2k = 0.0_wp
   ! !$omp end single
   ! !$omp do reduction(+:e_l2s,e_l2k,e_pl2s,e_pl2k)
   side_loop: do is=1,grid%ns
     xgs = affmap(grid%s(is),base%xigs)
     wgs = (grid%s(is)%a/base%me%voldm1)*base%wgs
     ! compute the diameters used in the mesh dependent norms
     hs = diameter( grid%s(is) )
     hk1 = diameter( grid%s(is)%e(1)%p )
     if(grid%s(is)%ie(2).gt.0) then
       hk2 = diameter( grid%s(is)%e(2)%p )
     else
       hk2 = 0.0_wp
     endif
     
     ! analytic solution
     lams = lam(xgs)
     ! L2 projection of the analytic solution
     do i=1,base%nk
       lamk(i) = sum( base%wgs*e_norm(i,:)*lams )
     enddo
     plams = matmul(lamk,e_norm)

     ! numerical solution
     lamhk = lamh((is-1)*base%nk+1 : is*base%nk)
     lamhs = matmul(lamhk,base%e)

     ! errors
     l2s  = sum(wgs*(lams-lamhs)**2)
     pl2s = sum(wgs*(plams-lamhs)**2)
     e_l2s  = e_l2s  + hs        * l2s
     e_l2k  = e_l2k  + (hk1+hk2) * l2s
     e_pl2s = e_pl2s + hs        * pl2s
     e_pl2k = e_pl2k + (hk1+hk2) * pl2s
   enddo side_loop
   ! !$omp end do nowait

   deallocate(  xgs,          &
     wgs, lams, plams, lamhs, &
     lamk, lamhk )
   ! !$omp end parallel
   deallocate( e_norm )

   if((e_l2s .lt.0.0_wp).or.(e_l2k.lt.0.0_wp).or.(e_pl2s.lt.0.0_wp).or. &
      (e_pl2k.lt.0.0_wp)) then
     write(message(1),'(a,e10.3,a,e10.3,a,e10.3,a,e10.3,a)') &
       'One error is negative: e_l2s**2 = ',e_l2s, &
       '; e_l2k**2 = ' ,e_l2k,  &
       '; e_pl2s**2 = ',e_pl2s, &
       '; e_pl2k**2 = ',e_pl2k,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_l2s  = max(e_l2s ,0.0_wp)
     e_l2k  = max(e_l2k ,0.0_wp)
     e_pl2s = max(e_pl2s,0.0_wp)
     e_pl2k = max(e_pl2k,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif
   e_l2s  = sqrt(e_l2s)
   e_l2k  = sqrt(e_l2k)
   e_pl2s = sqrt(e_pl2s)
   e_pl2k = sqrt(e_pl2k)
 
 end subroutine hybrid_err
 
!-----------------------------------------------------------------------

 !> Error computation for the primal variable.
 !!
 !! The following error is computed:
 !! <ul>
 !!  <li>\c e_l2
 !! \f$=\left( \sum_K \|u - u_h\|^2 \right)^{\frac{1}{2}} \f$
 !! </ul>
 !! where \f$K\f$ are the mesh elements.
 !!
 !! \note For consistency with the other subroutines in this module
 !! the format of \c uuuh is the same as for a discontinuous method,
 !! with element based degrees of freedom.
 !<
 subroutine primal_err(grid,base_in,uuuh,uuu,extra_deg, &
                       e_l2)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base_in
  real(wp), intent(in) :: uuuh(:)
  procedure(i_coeff_lam) :: uuu
  integer, intent(in) :: extra_deg
  real(wp), intent(out) :: e_l2

  integer :: ie
  real(wp) :: l2
  real(wp), allocatable :: xg(:,:), wg(:), uuue(:), uuuhk(:), uuuhe(:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'primal_err'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , p_s=base_in%p_s , &
                    deg=base_in%deg+extra_deg )
   allocate( xg(grid%m,base%m), wg(base%m), uuue(base%m), &
             uuuhk(base%pk), uuuhe(base%m) )

   e_l2  = 0.0_wp
   elem_loop: do ie=1,grid%ne
     xg = affmap(grid%e(ie),base%xig)
     wg = grid%e(ie)%det_b * base%wg

     ! analytic solution
     uuue = uuu(xg)

     ! numerical solution
     uuuhk = uuuh((ie-1)*base%pk+1 : ie*base%pk)
     uuuhe = matmul(uuuhk,base%p)

     ! errors
     l2  = sum(wg*(uuue-uuuhe)**2)
     e_l2 = e_l2 + l2

   enddo elem_loop
   deallocate( xg, wg, uuue, &
               uuuhk, uuuhe )

   if(e_l2.lt.0.0_wp) then
     write(message(1),'(a,e10.3,a)') &
       'One error is negative: e_l2**2 = ',e_l2,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_l2 = max(e_l2,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif

   e_l2  = sqrt(e_l2)

 end subroutine primal_err

!-----------------------------------------------------------------------

 !> Analogous to \c primal_err, but the numerical solution is given by
 !! the product of a polynomial function and an exponential.
 !<
 subroutine primal_err_sef(grid,base_in,uuuh,xi,xie,uuu,extra_deg, &
                           e_l2)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base_in
  real(wp), intent(in) :: uuuh(:)
  procedure(i_coeff_xiadv) :: xi
  real(wp), intent(in) :: xie(:)
  procedure(i_coeff_lam) :: uuu
  integer, intent(in) :: extra_deg
  real(wp), intent(out) :: e_l2

  integer :: ie
  real(wp) :: l2
  real(wp), allocatable :: xg(:,:), wg(:), uuue(:), uuuhk(:), uuuhe(:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'primal_err_sef'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , p_s=base_in%p_s , &
                    deg=base_in%deg+extra_deg )
   allocate( xg(grid%m,base%m), wg(base%m), uuue(base%m), &
             uuuhk(base%pk), uuuhe(base%m) )

   e_l2  = 0.0_wp
   elem_loop: do ie=1,grid%ne
     xg = affmap(grid%e(ie),base%xig)
     wg = grid%e(ie)%det_b * base%wg

     ! analytic solution
     uuue = uuu(xg)

     ! numerical solution
     uuuhk = uuuh((ie-1)*base%pk+1 : ie*base%pk)
     uuuhe = matmul(uuuhk,base%p) * exp(xi(xg)-xie(ie))

     ! errors
     l2  = sum(wg*(uuue-uuuhe)**2)
     e_l2 = e_l2 + l2

   enddo elem_loop
   deallocate( xg, wg, uuue, &
             uuuhk, uuuhe )

   if(e_l2.lt.0.0_wp) then
     write(message(1),'(a,e10.3,a)') &
       'One error is negative: e_l2**2 = ',e_l2,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_l2   = max(e_l2,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif

   e_l2  = sqrt(e_l2)

 end subroutine primal_err_sef

!-----------------------------------------------------------------------

 !> Error computation for the dual variable.
 !!
 !! The following errors are computed:
 !! <ul>
 !!  <li>\c e_l2
 !! \f$=\left( \sum_K \|\underline{q} - \underline{q}_h\|^2
 !!            \right)^{\frac{1}{2}} \f$
 !!  <li>\c e_hdiv
 !! \f$=\left( \sum_K \|\underline{q} - \underline{q}_h\|^2 + 
 !!            \sum_K \|\nabla\cdot\underline{q} -
 !!                     \nabla\cdot\underline{q}_h \|^2
 !!            \right)^{\frac{1}{2}} \f$
 !! </ul>
 !! where \f$K\f$ are the mesh elements.
 !!
 !! \note For consistency with the other subroutines in this module
 !! the format of \c qqqh is the same as for a discontinuous method,
 !! with element based degrees of freedom.
 !!
 !! \note The \f$H(div)\f$ is computed summing the element
 !! contributions: clearly, this only makes sense if
 !! \f$\underline{q}_h\f$ is globally \f$H(div)\f$.
 !!
 !! \bug The interface of \c qqq should be defined using the abstract
 !! interfaces given in \c mod_testcases.
 !<
 subroutine dual_err(grid,base_in,qqqh,qqq,div_qqq,extra_deg, &
                     e_l2,e_hdiv)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base_in
  real(wp), intent(in) :: qqqh(:)
  procedure(i_coeff_q) :: qqq
  procedure(i_coeff_divq) :: div_qqq
  integer, intent(in) :: extra_deg
  real(wp), intent(out) :: e_l2, e_hdiv

  integer :: ie, i
  real(wp) :: l2q, l2dq
  real(wp), allocatable :: xg(:,:), wg(:), &
    qqqe(:,:), dqqqe(:), qqqhk(:), qqqhe(:,:), dqqqhe(:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'dual_err'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , o_s=base_in%o_s , &
                    deg=base_in%deg+extra_deg )
   allocate( xg(grid%m,base%m),    wg(base%m),     &
           qqqe(grid%m,base%m), dqqqe(base%m),     &
             qqqhk(base%mk), qqqhe(grid%m,base%m), &
             dqqqhe(base%m) )

   e_l2   = 0.0_wp
   e_hdiv = 0.0_wp
   elem_loop: do ie=1,grid%ne
     xg = affmap(grid%e(ie),base%xig)
     wg = grid%e(ie)%det_b * base%wg

     ! analytic solution
      qqqe = qqq(xg)
     dqqqe = div_qqq(xg)

     ! numerical solution
     qqqhk = qqqh((ie-1)*base%mk+1 : ie*base%mk)
     do i=1,grid%d
       qqqhe(i,:) = matmul(qqqhk,base%o(i,:,:))
     enddo
     ! Piola transformation
     qqqhe = (1.0_wp/grid%e(ie)%det_b) * &
               matmul(grid%e(ie)%b,qqqhe(1:grid%d,:))
     ! divergence
     dqqqhe = (1.0_wp/grid%e(ie)%det_b) * matmul(qqqhk,base%divo)

     ! errors
     l2q  = sum(wg*sum((qqqe-qqqhe)**2,1))
     l2dq = sum(wg*(dqqqe-dqqqhe)**2)
     e_l2 = e_l2 + l2q
     e_hdiv = e_hdiv + l2q + l2dq

   enddo elem_loop
   deallocate( xg, wg,    &
             qqqe, dqqqe,  &
             qqqhk, qqqhe, dqqqhe )

   if((e_l2.lt.0.0_wp).or.(e_hdiv.lt.0.0_wp)) then
     write(message(1),'(a,e10.3,a,e10.3,a)') &
       'One error is negative: e_l2**2 = ',e_l2, &
       '; e_hdiv**2 = ',e_hdiv,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_l2   = max(e_l2,0.0_wp)
     e_hdiv = max(e_hdiv,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif

   e_l2   = sqrt(e_l2)
   e_hdiv = sqrt(e_hdiv)

 end subroutine dual_err

!-----------------------------------------------------------------------
 
 !> This subroutine computes the \f$L^2\f$ error for the total flux
 !! for a continuous primal formulation:
 !! <ul>
 !!  <li>\c e_flux_l2
 !! \f$=\left( \sum_K \|q - q_h\|^2 \right)^{\frac{1}{2}} \f$
 !! </ul>
 !! where \f$K\f$ are the mesh elements and the total flux \f$q_h\f$
 !! is computed as \f$-\mu\nabla u_h + au_h\f$.
 !!
 !! \note For consistency with the other subroutines in this module
 !! the format of uuuh is the same as for a discontinuous method, with
 !! element based degrees of freedom.
 !<
 subroutine primal_flux_err(grid,base_in,uuuh,qqq,coeff_diff,coeff_adv,&
   extra_deg, e_flux_l2 )
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base_in
  real(wp), intent(in) :: uuuh(:)
  procedure(i_coeff_q) :: qqq
  procedure(i_coeff_diff) :: coeff_diff
  procedure(i_coeff_adv ) :: coeff_adv
  integer, intent(in) :: extra_deg
  real(wp), intent(out) :: e_flux_l2
 
  integer :: ie, i, l
  real(wp) :: l2q
  real(wp), allocatable :: xg(:,:), wg(:), qqqe(:,:), mu(:,:,:),a(:,:),&
    uuuhk(:), uuuhe(:), grad_uuuhe(:,:), qqqhe(:,:)
  type(t_base) :: base ! basis with higher accuracy
  character(len=100) :: message(3)
  character(len=*), parameter :: &
    this_sub_name = 'primal_flux_err'

   ! increase the order of the quadrature formula
   base = new_base( base_in%me , p_s=base_in%p_s , &
                    deg=base_in%deg+extra_deg )
   allocate( xg(grid%m,base%m), wg(base%m),               &
     qqqe(grid%m,base%m), mu(grid%m,grid%m,base%m), &
     a(grid%m,base%m), uuuhk(base%pk), uuuhe(base%m),     &
     grad_uuuhe(grid%m,base%m), qqqhe(grid%m,base%m))

   e_flux_l2 = 0.0_wp
   elem_loop: do ie=1,grid%ne
     xg = affmap(grid%e(ie),base%xig)
     wg = grid%e(ie)%det_b * base%wg

     ! analytic solution
     qqqe = qqq(xg)

     ! problem coefficients
     mu = coeff_diff(xg)
     a  = coeff_adv(xg)

     ! numerical solution
     uuuhk = uuuh((ie-1)*base%pk+1 : ie*base%pk)
     ! scalar solution
     uuuhe = matmul(uuuhk,base%p)
     ! gradient of the scalar solution
     do i=1,grid%d
       grad_uuuhe(i,:) = matmul(uuuhk,base%gradp(i,:,:))
     enddo
     grad_uuuhe = matmul( transpose(grid%e(ie)%bi), &
                          grad_uuuhe(1:grid%d,:) )
     do l=1,base%m
       qqqhe(:,l) = -matmul(mu(:,:,l),grad_uuuhe(:,l)) &
                    + a(:,l)*uuuhe(l)
     enddo

     ! errors
     l2q  = sum(wg*sum((qqqe-qqqhe)**2,1))
     e_flux_l2 = e_flux_l2 + l2q

   enddo elem_loop
   deallocate( xg, wg, qqqe,          &
             mu, a,                    &
             uuuhk, uuuhe, grad_uuuhe, &
             qqqhe)

   if(e_flux_l2.lt.0.0_wp) then
     write(message(1),'(a,e10.3,a)') &
       'One error is negative: e_flux_l2**2 = ',e_flux_l2,';'
     write(message(2),'(a)') &
       '  A possible reason is the use of a quadrature formula with'
     write(message(3),'(a)') &
       '  negative weights. Error set to zero.'
     e_flux_l2 = max(e_flux_l2,0.0_wp)
     call warning(this_sub_name,this_mod_name,message)
   endif

   e_flux_l2 = sqrt(e_flux_l2)
 
 end subroutine primal_flux_err
 
!-----------------------------------------------------------------------

end module mod_error_norms

